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ABSTRACT 

SCORCH (Simulations and Constructions of the Reionization of Cosmic Hydrogen) is a new project to 
study the Epoch of Reionization (EoR). In this first paper, we probe the connection between observed 
high-redshift galaxies and simulated dark matter halos to better understand the primary source of 
ionizing radiation. High-resolution N-body simulations are run to quantify the abundance of dark 
matter halos as a function of mass M, accretion rate M, and redshift z. A new fit for the halo 
mass function dn/dM is ~ 20% more accurate at the high-mass end. A novel approach is used to 
fit the halo accretion rate function dn/dM in terms of the halo mass function. Abundance matching 
against the observed galaxy luminosity function is used to estimate the luminosity-mass relation and 
the luminosity-accretion-rate relation. The inferred star formation efficiency is not monotonic with 
M nor M, but reaches a maximum value at a characteristic mass ^ 2 x 10^^ Mq and a characteristic 
accretion rate ^ 6 x 10^ M0yr“^ at z ~ 6. We find a universal EoR luminosity-accretion-rate relation 
and construct a fiducial model for the galaxy luminosity function. The Schechter parameters evolve 
such that (/>* decreases, M* is fainter, and a is steeper at higher redshifts. We forecast for the upcoming 
James Webb Space Telescope and show that with apparent magnitude limit toab ~ 31 (32), it can 
observe >11 (24) unlensed galaxies per square degree per unit redshift at least down to M* at z < 13 
(14). 

Keywords: cosmology: theory - dark ages, reionization, first stars - galaxies: high-redshift - large-scale 
structure of universe - methods: numerical 


1. INTRODUCTION 

Cosmic reionization is a frontier topic in cosmology 
with plenty of scientific richness for theoretical and obser¬ 
vational explorations. The epoch of reionization (EoR) is 
uniquely marked by the emergence of the first luminous 
sources: stars, galaxies, and quasars in the first billion 
years. It is possibly the only time that luminous sources 
directly and dramatically alter the state of the entire 
universe, converting the cold and neutral intergalactic 
medium (IGM) into a warm and highly ionized one. 
Studying the EoR will reveal how the first generation 
of stars, galaxies, and black holes formed and evolved. 
It can also provide complementary constraints on cos¬ 
mological parameters similar to stu dies of the cosmic 
microw ave background (CMB). See iLoeb &: Furlanett^ 
(|2013ll for a recent review. 

Current observations suggest that reionization was al¬ 
ready in significant progress by z ^ 9 and must have 
ended by z ^ 6. The photoionization and photo¬ 
heating of hydrogen leaves many observable imprints. 
Neutral hydrogen can be detected t hrough 21cm radi¬ 
ation from radio observat ions (e.g. iBowman fc RogersI 
I2OIOI : iParsons et al.1 1201^. in abs o rption as th e Ly¬ 
man alpha forest le.g. IPan et al.l 120021 1200011 . and 
also through Lyman alpha scatter i ng of high-redshif t 
galaxies (e.g. iKashikawa et al.l l2006t lOuchi et al.ll2009ll . 
Free electrons will scatter CM B photons and produce 
patchy Thomson scattering (e.g . iDvorkin fc Smithll2009l : 
iSii et al.l I2011I iNataraian et al.l 1201311 . kinetic Sunyaey- 
Zel’do yich temperature anis o tropy o n arcminute scales 
(e.g. iSiinyaev fc Zeldovichl I1970I iZahn et ^ 120121 : 


iSieyers et al.l 1201311. and oolMizat i on anisotropy at low 
multipoles fe.g. iK ogiit et ^120031 iHinshaw et al.ll2?)T^ 
iPlanck Collaboration et al.l 120151 1 Heating of the IGM 
affec ts the thermal broadening of the Lyman alpha forest 
fe.g. iTheuns et ani2002t iGen et aDl2009t iLidz fc Malloyl 

1201411 . Ongoing and upcoming obseryations will pro- 
yide multi-wavelength constraints on the uncertain tim¬ 
ing and duration of the EoR. 

Current theory suggest that large-scale, overdense re¬ 
gions near radiation sources are generally reionized ear¬ 
lier than large-scale, underdense regions far from sources. 
Thre e main approaches a re used to model the EoR 
(e.g. iTrac fc GnedinI [20111 1. Gosmological simulations 
combining N-body, hydro, and radiative transfer (RT) 
algorithms are the most accurate, but also most ex¬ 
pensive approach to solve the coupled ev olution of the 
dark matter, baryons, and radiation fe.g iGnedinI [20001 : 
iTrac et ^120081 : lGnedinir20l4 iNorman et al.ll2015ll . RT 

post-processed on saved N-body or hydro simulated 
data can be more cost effective while capturing im- 
portant features of nonlinear structure formation (e.g . 

Ciardi et al.ll2003Hlliev et al.ll200(THMcQuinn et al.ll2007l : 

Finlator et al.l 120091 : lAubert fc Tevssied 1201011 . Semi- 
analytical/numerical methods provide an approximate 
and efficient approach to exp lore the vast parame¬ 
ter space o f reio n ization (e.g. iFiirlanetto et al.l 12004 
Zahn et al.l 120071: iSantos et al.lT2010l : iMesinger et al.l 
201lHBattaglia et al.ll2013bl l. Recently, there is renewed 
emphasis on using radiation-hydrodynamic simulations 
to model the complex astrophysics of sources and sinks 
and using semi-numerical methods to make theoreti¬ 
cal predictions and mock observations. These two ap- 
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proaches in tandem provide our best option in making 
forward progress in synergy with observations. 

SCORCH (Simulations and Constructions of the 
Reionization of Cosmic Hydrogen) is a new project to 
study the EoR and provide useful theoretical tools and 
predictions to facilitate more accurate and efficient com¬ 
parison between observations and theory. To build a 
robust theoretical framework, we need to further in¬ 
vestigate how the distribution and properties of radi¬ 
ation sources and sinks affect the complex reioniza¬ 
tion process. Th e RadHydro code (|Trac fc PenI 12004 
iTrac fc CenI 120071 1 will be used to produce N-body -I- 
hydro -|- RT simulations with subgrid physics modeled 
using the latest observations and simulations. Mock 
observations will be constructed by mapping higher- 
resolution, smaller-volume radiation-hydrodynamic sim¬ 
ulations onto l ower-resolution, larger-volume N-body 
simulations (e.g iBattaglia et al.ll2013tj: iNataraian et al.l 


l2013t IBattaglia et al.ll2013aHLa Plante et al.ll2014 K 
Currently, there are ~ 1500 galaxy candidates ob¬ 
served in the redshift range 6 < z < 10 from Hubble 
Spac e Telescope (HST) observ a tions (e.g. iGrogin et al l 
2011|_jKoe cemoer^^^ 1 120111: IWindhorst et al.1 120111 


Ellis ^et aiT 201.1l : llllingworth et al.l 120131 : iSchmidt et alJ 

2014f) . The galaxy luminosity function is well fit 


by a Schechter function and the redshift depen- 
dence of the parameters have been quantified (e.g. 
Bouwens et al.l [20 H: lOes ch et al.l 12014 iBouwens et al.l 


2015at iFinkelstein et al.ll2015al l. However, it is unclear 


what is the physical origin for the evolution and how 
is it connected to the growth of structure, particularly 
that of dark matter halos. The abundance of dark mat¬ 
ter halos ar e accurately quanti f ied using N-body simu- 
lations (e.g. iJenkins et al.l 120011 : 1 Tinker e t al.ll200j ^. biR 
previous work rel ated to the EoR, te.g lReed et al.ll2007l : 
iLnkic et'm 120071 1 have not extensively studied atomic 
cooling halos (T > 10^ K, M > 10® Mq) capable of host¬ 
ing high-redshift galaxies. Similarly, the growth of dark 


matter halos have also been measured using N-body sim¬ 
ulations fe.glWechsler et al.ll2002t iFakhouri fc Mall2008l : 
iCorrea et al.ll2015bh . but generally for higher mass and 
at lower redshift. 

In Paper I, we quantify the connection between ob¬ 
served high-redshift galaxies and simulated dark mat¬ 
ter halos in order to better understand the abundance 
and evolution of the primary source of ionizing radiation. 
Section [5] describes the N-body code and halo finder used 
to simulate and locate dark matter halos. Sections [3] and 
in describe how the abundance of dark matter halos as a 
function of mass M, accretion rate M, and redshift z are 
quantified. In Sections [5] and El we use the abundance 
matching technique to estimate the luminosity-mass rela¬ 
tion and the luminosity-mass-accretion-rate relation and 
infer the star formation efficiency. In Sections 0 and [H 
we make predictions for the high-redshift galaxy lumi¬ 
nosity function and forecast galaxy counts for the up¬ 
coming James Webb Space Telescope (JWST). Section |9] 
discusses the implications of our results on cosmic reion¬ 
ization and Section ITUI concludes with a summary of ma¬ 
jor results. We adopt the concordance cosmological pa¬ 
rameters: = 0.27, Ha = 0.73, Hb = 0.045, h — 0.7, 

(78 = 0.8, and rig = 0.96. 


2. N-BODY SIMULATIONS 

N-body simulations are run using a new 
particle-particle-partic l e-mes h code (P®M; e.g. 
iHocknev fc EastwoodI 119811) containing updated al¬ 
gorithms. The gravitational potential is decomposed 
into short and l ong-ra.nge co mponents using a gaussian 
kernel following iBaglal (I2002D . The long-range potential 
is computed using a particle-mesh (PM) algorithm 
where Poisson’s equation is efficiently solved using Fast 
Fourier Transforms (FFTs). Particles are assigned to a 
mesh using the triangular-shaped-cells (TSC) scheme 
to reduce mesh anisotropy. The short-range force is 
computed for particle-particle (PP) interactions using 
direct summation, which is only moderately expensive 
at high redshifts when the large-scale structure is not 
highly clustered. Adaptive time ste pping is perform ed 
using the kick-drift-kick scheme from iSoringell (|2005D . 

Dark matter halos are found using a new hybrid finder 
that has two ma jor components. F irst, a friends-of- 
friends (FoF; e.g. I.Ienkins et al.ll2nnil) algorithm is used 
to locate the peaks of halo candidates. A linking length 
h = 0.08 times the mean interparticle separation is cho¬ 
sen to pick o ut only the inner regio ns and prevent over¬ 
merging (e.g. iCohn fc White! l2008l) . For each peak, the 
particle with the minimum value of the density-weighted 
gravitational potential pcf) is chosen as the center. _ 

Sec ond, a spherical overdensity (SO; e.g. lLacev fc Colei 
11994) algorithm is used to measure halo properties. A 
halo with mass Ma within radius Ra has an average 
density that is factor of Ahaio times the average matter 
density pm such that 

Ma = ^7r(AhaloyOm).RA- (1) 

Throughout the paper, M and R refer to M 200 and i? 200 ) 
respectively. Halos are allowed to overlap, but if two can¬ 
didates have centers within the larger halo’s radius, one 
candidate is chosen as the primary halo and the other is 
considered a subhalo. The primary halo is the one that 
satisfies more of the following criteria: larger masses and 
more negative potentials within i ?200 and i? 5 oo- Only 
halos with at least 400 particles (l/-\/400 = 0.05) are 
counted to avoid resolution effects and ensure complete¬ 
ness near the low-mass end of each simulation. 

The halo finder is run “on the fly” every 20 million cos¬ 
mic years. This time interval is comparable to relevant 
timescales such as the halo dynamical time and the mass 
accretion timescale at high redshifts. In constructing the 
halo merger trees, every particle is tracked rather than 
just the most-bound particle in identified halos. We use 
both forward and backward linking to carry out a two- 
step process in progenitor-descendent matching. First, 
for a given progenitor of mass M 2 at a given time t 2 , its 
descendent of mass Mi from a previous time ti is the one 
that contributes the largest weight, which is defined to be 
the sum of the gravitational binding energy from parti¬ 
cles in the progenitor that once belonged to the descen¬ 
dent. Second, the selected descendent may contribute 
particles to other halos, but the weighted contribution 
must be less than what it gave to the main progenitor. 
The mass accretion rate is calculated as 


t2 — ti 


(2) 
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Table 1 

N-body simulation parameters 


Name 

L 

(Mpc h~^) 

nip 

{Meh 

-1) 

e 

(kpc h~^) 

^init 

-^halo,min 

(Mq/i-I) 

Mace, 

{Meh 

min 

NB_L10.N2048 

10 

8.72 

X 

10® 

0.31 

350 

3.49 

X 

10® 

2.18 

X 

lO’’ 

NB_L15_N2048 

15 

2.94 

X 

10^ 

0.46 

350 

1.18 

X 

107 

7.36 

X 

107 

NB_L25.N2048 

25 

1.36 

X 

10® 

0.76 

300 

5.45 

X 

10" 

3.41 

X 

10® 

NB_L35.N2048 

35 

3.74 

X 

10® 

1.07 

300 

1.50 

X 

10® 

9.35 

X 

10® 

NB_L50.N2048 

50 

1.09 

X 

10® 

1.53 

300 

4.36 

X 

10® 

2.73 

X 

10® 

NB_L75.N2048 

75 

3.68 

X 

10® 

2.29 

300 

1.47 

X 

lO'^ 

9.20 

X 

10® 

NB_L100.N2048 

100 

8.72 

X 

10® 

3.05 

250 

3.49 

X 

10® 

2.18 

X 

lOi® 

NB_L150_N2048 

150 

2.94 

X 

107 

4.58 

250 

1.18 

X 

IQio 

7.36 

X 

lOlo 

NB_L200.N2048 

200 

6.98 

X 

10^ 

6.10 

200 

2.79 

X 

IQi® 

1.74 

X 

loll 

NB_L300.N2048 

300 

2.36 

X 

10® 

9.16 

200 

9.42 

X 

IQlO 

5.89 

X 

loll 

NB_L400.N2048 

400 

5.58 

X 

10® 

12.2 

150 

2.23 

X 

IQll 

1.40 

X 

1012 


Note. — Two realizations of each P^M N-body simulation containing 2048^ dark matter particles are run using 2LPT initial conditions. 
Dark matter halos are found with a hybrid algorithm and have an average density phalo = 200pm- For calculating the halo mass and 
accretion rate functions, the minimum masses correspond to 400 and 2500 particles, respectively. 


Since it is more difficult to robustly measure accretion 
rate than mass, only halos with at least 2500 particles 
(1/V2500 = 0.02) are used. 

Table [1] lists the N-body simulations used to quan¬ 
tify the abundance of halos as a function of mass, ac¬ 
cretion rate, and redshift. Box sizes in the range 10 < 
L/(Mpc/i“^) < 400 are chosen to focus on atomic cool¬ 
ing halos and two realizations of each box size are run to 
reduce sample variance. The nine largest box-size simula¬ 
tions are necessary to reach convergence at the lowest M 
of interest, while the two smallest box-size simulations 
are for convergence at the lowest M. Each simulation 
contains 2048^ dark matter particles and has a particle 
mass resolution 8.72 x 10®(L/100)^ The grav¬ 

itational softening length is set to 1/16*^'' of the mean 
interparticle spacing. Initial conditions are generated 
using 2nd-order la grangian perturbation theory (2LPT; 
lScoccim arro||19 98f) with a lin ear transfer function from 
CAMB (jb^wis fc Bridlel[2?)?)^ . The starting redshift Zinit 
is chosen such that the perturbations, (5rms ~ 0.01, are 
still linear. 

3. HALO MASS FUNCTION 


(|2006ll : [Tinker et al.l (I2008D . the barrier-crossing distri¬ 
bution function is parametrized as 


/(o-) = A 



e 


-c/ 


2 


cr 


( 5 ) 


where A sets the overall amplitude, a and b set the slope 
and amplitude of the low-mass power law, and c sets the 
exponential decrease scale. 

Missing large-scale power due to f inite box sizes ar e 
corrected using a similar procedure to lReed et al.l (1200711 . 
An effective variance is calculated as 

als{M,z) = ^ \Sik,z)f\Wik,M)\\ (6) 

k 


where i5(k, z) is the Fourier transform of the linear over¬ 
density field extrapolated to redshift z. An effective mass 
Meff is calculated by finding the mass in Equation 0] 
which would yield a variance equal to The effective 
quantities and tJeff are used instead of their normal 
counterparts to calculate the halo mass function. 

We combine all of the simulations and bin the halos to 
calculate the discrete halo mass function. 


The halo mass function in differential form is defined 
as the comoving number density n(> M,z) of halos of 
mass M per unit mass dM. In extended Press-Schechte r 
theory (EPS; e.g. lBond et al.l[T9MI : iLacev fc Colelll993l l. 
it is generally expressed as 


dn 

dM 


fi<^) 


po din a ^ 
M dM 


( 3 ) 


where po is the comoving average cosmic density, a is 
the rms fluctuation of the smoothed linear density field, 
and /(cr) is the cr-weighted distribution of random-walk 
barrier-crossings. For a continuous density field, the vari¬ 
ance of the density fluctuations is calculated as 



P{k,z)\W{k,M)\^k^dk, 


( 4 ) 


where P{k, z) is the linear power spectrum extrapo¬ 
lated to redshift z and W{k,M) is the Fourier trans¬ 
form of the spherical tophat window function of ra¬ 
dius R = [M/(4/37rpo)]^^^- Following iWarren et al.l 


=y . (7) 

AVAM ^ ViAM ^ ’ 

i 

Each binned halo is given a weight w = M/M^ff such 
that total mass in any given mass bin is conserved. Since 
simulations with different box sizes have different mini¬ 
mum masses, the effective volumes V differ for the vari¬ 
ous mass bins. 

Figure 0] shows the differential mass functions and 
barrier-crossing distribution functions at z = 6, 8, and 
10. The mass functions have a generic shape with a low- 
mass power-law and a high-mass decaying tail. The func¬ 
tion /(cr) shows no redshift dependence after correcting 
for finite box-size effects. Since our box sizes are reason¬ 
ably large for the mass and redshift ranges considered, 
the corrections are typically small, only ~ 10% at z = 10 
and ^ 1% at z = 6 for moderate masses. 

The middle panels of Figure [T] show a comp arison with 
the best-ht results from [Tinker et al.l (1200811 , which are 
calibrated for the mass range 10^^ < M/{MQh~^) < 
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Figure 1. Left: The differential halo mass function at 2 = 6 (blue squares), 8 (green circles), and 10 (red triangles) from a series of N-body 
simulations (top). Comoving number densities are in units of Mpc~^h^ a nd ha lo masses are in units of The simulation results 

are in good agreement overall with the best-fit results from [Tinker et HI II2008I) . but differences of 20% are seen at the high-mass end 
(middle). For our new fit, the residuals are typically < 5% across the mass and redshift ranges of interest (bottom). Right: The barrier¬ 
crossing distribution function /((j) has a universal form for the scales and redshifts considered. The variable <j~^ increases monotonically 
with mass. Poisson error bars are shown if they are larger than the data point symbols. 


10^® and redshift range 0 < z < 2.5. As suggested in 
their paper and through personal communication, we use 
a maximum redshift z = 2.5 in their equations for the 
redshift evolution of the fitting parameters to obtain A = 
0.156, a = 1.36, b = 2.54, and c = 1.19. The agreement 
is remarkably good overall considering the extrapolation 
in both mass and redshift. It is better at lower masses 
and at lower redshifts, but differences of ~ 20% are seen 
at the high-mass end where bright high-redshift galaxies 
are expected to reside. 

Since f{a) has a universal form for the scales and red¬ 
shifts considered, we choose to fit the z = 6 results only 
because of the larger range and higher signal-to-noise. 
Furthermore, the parameters a and b are kept unchanged 
because the simulated halo catalogs do not sample the 
high-(T (low-mass) power law portion of the function. 
The best-fit barrier-crossing distribution function is 


/(a) =0.150 



( 8 ) 


The bottom panels of Figured] show that the residuals for 
the new fit are < 5% across the mass and redshift ranges 
of interest. The larger deficits at the highest mass and 
particularly at higher redshifts arise because we cannot 
perfectly correct for the missing large-scale power in fi¬ 
nite simulation boxes. In the exponential tail of the halo 
mass function, a small mass change can lead to a rela¬ 
tively large change in the number density. While we have 
started the simulations at high enough initial redshifts 
such that the second-order displacement corrections are 
themselves small, starting at even high redshifts may re¬ 
duce the mass deficits. 

The halo mass function can be calculated for any cos¬ 


mological model given the linear power spectrum. For 
our set of cosmological parameters, the tophat variance 
can be fit to < 2% accuracy for lO’’ < M/{MQh~^) < 
10^3 with 


cr i(z) = D i(z) 


0.11-f 0.077 


M 


108 


0.18 


(9) 

where D{z) is the linear growth factor. Note that a 
double power law plus a constant term can accurately 
fit the entire mass range spanning mini-halos (M ~ 
10® to massive clusters {M ~ 10^® MQh~^). 


4. MASS ACCRETION RATES 

The halo accretion rate function in differential form is 
defined as the comoving number density n{> M, z) of 
halos with accretion rate M per unit accretion rate dM. 
In principle, its functional form in terms of accretion rate 
and redshift can be motivated with EPS theory. Alter¬ 
natively, we take a novel approach and express it as 

dn dn dM 
dM ~ dM dM' 

where dn/dM is the differential mass function (Equa¬ 
tion [3|). The advantage of using this functional form is 
that the redshift dependence of the halo mass function 
is already well known. To use Equation (TO] we need to 
understand how mass and accretion rate are related and 
be able to calculate the derivative dM / dM. 

Eigure [5] shows how the mass accretion rate depends 
on mass at z = 6, 8, and 10. The average accretion rates 
for the mass range 10® < M < 10^® and redshift range 
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Figure 2. Left: The average mass accretion rates at 2 = 6 (blue squares), 8 (green circles), and 10 (red triangles) are accurately fit by 
a (M) oc M^-^®(1 + 2 )^-^ relation (top). The SO halo rates are generally higher than the FoF halo rates from IFakhouri et al.l II201QI) . but 
there is better agreement in the overlap mass range M > 10^^ Mq (middle). For our new fit, the residuals are typically < 5% across the 
mass and redshift ranges of interest (bottom). Right: The variance and skewness also have power-law relations with mass, but with slightly 
shallower slopes and weaker redshift dependences. The statistical uncertainties in the binned mean values are shown if the error bars are 
larger than the data point symbols. 


6 < z < 10 are well fit by 


(M) 


0.21MQyr-i 


M 

108 Mq 



( 11 ) 


The residuals are typically < 5% across the mass and red¬ 
shift ranges of interest. The larger deficits at the highest 
masses and particularly at higher redshifts are related to 
those seen in the halo mass function in Figure [T] since the 
accretion rates are also affected by the missing large-scale 
power in finite simulation boxes. 

Our mass power-law slope of 1.06 is slight l y shal lower 
than the value of 1.1 from IFakhouri et al.l (I 2010 D and 
slightly steeper than the linear mass dependence from 
iCorrea et al.l ()2015al lbll. We find similar results for M > 
10 ^*^ Mq where our mass range overlaps with theirs, but 
our accretion rates become relatively larger (smaller) at 
lower mass because of the shallower (steeper) slope. We 
all find the same r edshif t power-law slope of 2.5 at high 
z. IFakhouri et al.l (|2010f) measure accretion rates from N- 
body simulations, but they use a FoF halo finder instead. 
Since the ratio of M 200 and MpoP is neither constant with 
mass nor redshift, it is not surprising to find differences 
between M 200 and MpoP- Their fit is done at low redshift 
z ^ 0 and it is not clear how well their best fit compar es 
to their own data for 6 < z < 10 . ICorrea et al.l (|2015al lbfl 
use EPS theory and N-body simulations to derive and 
measure the mass accretion history, respectively. They 
find a linear dependence on mass from their analytical 
work and assume this same mass dependence when fitting 
the numerical results. Thus, it is not clear if the latter 
would have had a slightly higher mass power-law slope. 
Nonetheless, there is overall good agreement between all 


of the results given the differences in details. 

The distribution of mass accretion rate at any given 
mass is positively skewed. The best-fit relations for the 
variance /r 2 = ((M — (M))^) and skewness ^3 = ((M — 
(M)) 8 ) are 


= 0 . 28 MQyr-i 

/iJ/^=O.42M0yr-i 


/ M + 


V10® Mq 
/ ^ X 1.0 

(108 M© 


1 + z 


0.6 


( 12 ) 

(13) 


Equations [T1][T3] can be used with, for example, an Edge- 
worth expansion to model the probability distribution 
function of accretion rates at a given mass and redshift. 
Mergers are responsible for the positive tail of the dis¬ 
tribution. Thus, the similarity in mass dependence and 

• 1 /o 

difference in redshift evolution between (M), , and 

1 /3 

/Tg reflect those between the average accretion rate and 
the average merger rate. See IFakhouri et al.l (|2010ll for 
recent work on halo merger rates. 

Figure |3] shows the abundance of dark matter halos as 
a function of mass accretion rate at z = 6 , 8 , and 10. The 
halo accretion rate function and the halo mass function 
are similar in shape, with each having a low-end power 
law and a high-end exponential decline. We can model 
the halo accretion rate function using Equation [10] and 
by combining the halo mass function with the best-fit 
mediating mass. 


M = 4.2 X 10®Mq 


M 

M0yr-i 


0.91 



(14) 


-2.4 
















































6 


Trac, Cen, & Mansfield 



1 10 ' 102 lo" 10 '' 

M [Me/yr] 


Figure 3. The differential halo accretion rate function at 2 : = 6 
(blue squares), 8 (green circles), and 10 (red triangles) from a series 
of N-body simulations (top). Comoving number densities are in 
units of Mpc~^ and mass accretion rates are in units of AfQyr”^. 
The halo mass function and the halo accretion rate function are 
similar in shape, with each having a low-end power law and a 
high-end exponential de clin e. T he l atter can be predicted from the 
former using Equations 1101 and 1141 The residuals are typically < 
10% for the mass and redshift ranges of interest (bottom). Poisson 
error bars are shown if they are larger than the data point symbols. 

The bottom panel of Figure[3]shows that the residuals for 
the best fit are typically < 10% for the mass and redshift 
ranges of interest. The larger deficits at the highest M 
and particularly at higher redshifts are related to those 
seen in the halo mass function in Figure [T] since the ac¬ 
cretion rates are also affected by the missing large-scale 
power in finite simulation boxes. Note that Equation 
M is not simply obtained from Equation IIII because the 
distribution of accretion rate at a given mass is skewed. 
However, they both have the convenient property of sep¬ 
arable power-law dependences on M, M, and z. 

The halo accretion rate function does not have a sim¬ 
ple redshift dependence according to Equation [101 For 
a given accretion rate, both the mediating mass and its 
derivative increase at lower redshifts. This is a conse¬ 
quence of accretion slowing down as seen in Figure ID 
In contrast, the halo mass function at a higher mass 
decreases in amplitude because more massive halos are 
rarer in hierarchical structure formation. We expect the 
redshift dependence in Equations ITTlfMl to hold for red¬ 
shifts z > 6 relevant to the EoR. However, our single 
redshift power law scaling is n ot appropriate at lower 
redshifts. iMcBride et al.l ( 200911 have shown that a more 
complex scaling is required to parametrize the entire 
mass accretion history. More work is required to exam¬ 
ine the validity and accuracy of extrapolating our results 
to lower redshifts. 


5. ABUNDANCE MATCHING 


Abundance matching is performed by equating the 
number density of galaxies to the number density of halos 
in differential form, 

dUgal du-halo dA I 

dLuy dX dLjjy ’ 
or in cumulative form. 


TuVi ■^) — 11'halo(^ A, 2^). (Hi) 


If X is taken to be the halo mass M, then we infer 
the luminosity-mass relation z). This procedure 

does not account for scatter in the mass-to-light ratio 
and the relation should be interpreted as an average lu¬ 
minosity for a given mass. If X is taken to be the mass 
accretion rate M, then we infer the luminosity-accretion 
rate relation Ltjy{M,z). This procedure accounts for 
the episodic nature of star formation and scatter in the 
mass-to-light ratio. Furthermore, it is more logical since 
the star formation rate should be more highly correlated 
with the halo growth rate rather than halo mass. 

The galaxy luminosity function is generally 
parametrized with a Schechter function in luminos¬ 
ity form. 


4>{L\jy) = (j) 

or in magnitude form, 




\ J 


(j){Muy) = (0.41nl0)(^* 




- cn-j-l 


2 ^q 0 . 4 (M,-Muv) 

X exp ['ioO- 4 (m*-Muv) 


(17) 


(18) 


where is an overall normalization, L* is a character¬ 
istic luminosity, M* is a characteristic magnitude, and a 
is the slope of the faint-end power law. The conversion 
between UV magnitude and luminosity is given by the 
standard AB relation. 


M\jy = —2.5 log 


L\jy 


4.487 X 10^° ergs“iHz“i 


(19) 


The cumulative galaxy luminosity function is obtained 
by integrating Equation [T7| analytically to get the in¬ 
complete gamma function (/)*L*F(1 -|- a, Luy/L^,). The 
cumulative halo mass function and mass accretion rate 
function do not have simple analytical forms, but are 
e asily obtained b y numer ical integration. 

iBouwens et al.l (l2015af) have made the most recent 
measurements of the high-redshift UV luminosity func¬ 
tions. For three redshift samples with (z) Ri 5.9, 
7.9, and 10.4, there are 867, 217, and 6 galaxy candi¬ 
dates observed, respectively. For reference, their best-fit 
Schechter parameters are: 


(z) « 5.9, 
(z) « 7.9, 
(z) 10.4, 

z ^ 10.4, 


(/)* = 5.0 X 
= 2.1 X 
= 8.0 X 
= 3.0 X 


10-^ M, 
10 -^ 
10-^ M* 
10-^ M* 


= -20.94, 
= -20.63, 
= -20.92, 
= -20.92, 


a = -1.87 
a = -2.02 
a = -2.27 
a = -2.27 


where </>* has units of comoving Mpc Note that the 
last entry above is not a best ht to the galaxy counts, but 
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Figure 4. Left: The UV magnitude as a function of halo mass for z ~ 6 (blue), 8 (green), and 10 (red) from abundance matching Luv 
and M. The shaded curves are illustrative of the la (68%) uncertainty in the currently-observed galaxy luminosity functions and are 
truncated at the star formation limit. Right: The UV magnitude as a function of mass accretion rate from abundance matching Luv 
M. The luminosity-accretion relation is more consistent with no redshift evolution than the luminosity-mass relation for the redshift range 
6 < ^ < 10. The luminosity monotonically increases with mass and accretion rate, but the effective power-law slopes suggest that the 
relations should be fit with, for example, triple power law functions (bottom). 


an extrapolation of the luminosity functions from lower 
redshifts. They allowed the parameters log(/>*, M*, and 
a to vary linearly with 2 and fitted the galaxy counts in 
the redshift range 4 < 2 < 8 in order to estimate the 
results for higher redshifts. Using the derived redshift 
dependence, they were able to fix M* and a and then 
fit for 0* at (z) « 10.4. We note that M* = —20.92 
at z ~ 10 is rather bright considering that the best-fit 
characteristic magnitude changed from -20.94 to -20.63 
for z Ri 6 to 8. Furthermore, the extrapolation from a 
different redshif t range 5 < z < 8 yields M* k. —20.22 
(|Bouwens et al.l [20151)1) . which is more in line with the 
best-fit values for z r; 6 — 8. 

To estimate upper and lower bounds on a given galaxy 
luminosity function, we vary the three Schechter parame¬ 
ters using their uncertainties, add the weighted Schechter 
functions to the M\jy-(j) grid, and calculate the 68% confi¬ 
dence level. For the weight, we use a trivariate gaussian 
likelihood that assumes uncorrelated errors. However, 
there is degeneracy between the parameters and their er¬ 
rors are correlated. To prevent overestimating the confi¬ 
dence region, we reduce the errors in the Schechter pa¬ 
rameters to 2/3 of their original values. This procedure 
allows us to match the Icr errors in the binned measure¬ 
ments of the galaxy luminosity function. A more rigorous 
statistical analysis would require using the full likelihood 
for the galaxy luminosity function fit, which is not avail¬ 
able. 

Figure |4] shows how the UV luminosity and magnitude 
depend on mass and accretion rate at z Ri 6, 8, and 10 
(i.e. 5.9, 7.9, and 10.4). The shaded curves are illustra¬ 
tive of the 1(7 (68%) uncertainty in the current galaxy lu¬ 
minosity function. Our error analysis is able to capture 


the trend of increasing uncertainty at higher redshifts. 
Neither the luminosity-mass relation nor the luminosity- 
accretion-rate relation is represented by a single power 
law, but appears to be composed of several power law 
portions. 

To qualitatively understand the shape of the 
luminosity-mass and luminosity-accretion-rate relations, 
we also plot their effective power-law slopes in the bottom 
panels of Figure 01 The cumulative galaxy and halo num¬ 
ber densities can be written in effective power form as 
ngai(L) oc nhaio(A7) oc and nhaio(-/U) oc 

M'l+'yeH _ luminosity-mass and luminosity-accretion- 
rate relations can be written as Luv oc and Luv oc 

From Equation [T51 the slopes are related through 


d In Luv 
“ d\nM 


dinrthaio / dlnugai 

dlnM / din Luv 


f + i^eff 

1 T CTeff 


( 20 ) 


d In Luv _ din rthaio / dlnugai __ 1 + 7eff . 21 ) 
dIn M dIn M ' din Luv 1 + oteS 


At low LuV) M, and M, where the galaxy luminosity 
function and halo mass and accretion rate functions are 
approximately power laws (aeff ~ ctj PeS ~ ~2, 7 eff Ri 
—2), the effective power-law slopes of the luminosity- 
mass and luminosity-accretion-rate relations are approxi¬ 
mately constant (5eff ~ Ceff ~ —1/(1 -I-q;). At intermediate 
scales, the halo mass and accretion rate functions deviate 
from their power law forms before the galaxy luminosity 
function does, both /3eff and 7 eff become more negative, 
and both deff and Cefi increase. As the galaxy luminosity 
function deviates from its power law form, aeff becomes 
more negative, and both des and decrease again. At 
the very bright and massive end, both the galaxy and 
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halo number densities are exponentially decreasing, but 
such that tteff < /3eff and < 7eff, resulting in Jeff and 
Ceff being approximately constant. The effective power- 
law slopes suggest that both relations can be fit with, for 
example, a triple power law. 


Luv = Lq 





c—b 

, ( 22 ) 


where Lq is an overall amplitude, a, 6, and c are three 
power-law slopes, and Xg,, Xb, and Xc are three charac¬ 
teristic scales. 

For the luminosity-mass relation, brighter galaxies are 
found in more massive halos at any given redshift as ex¬ 
pected from the abundance matching procedure. At any 
given mass, brighter galaxies are found at higher redshifts 
and this is a consequence of the relative evolution of the 
galaxy lumino sity and halo mass function s . Our results 
are similar to iKuhlen &: Faucher-Giguerel (120121 1. where 
they performed abundance matching using older fits to 
the galaxy luminosity function and halo mass function. 

For the luminosity-accretion-rate relation, brighter 
galaxies are found in halos with larger accretion rates 
at any given redshift by construction. The z ~ 6 — 8 re¬ 
sults are remarkably very similar for the entire accretion 
rate of interest. In fact, the luminosity-accretion-rate 
relation is consistent with no evolution for the redshift 
range 6 z < 10 when taking into account the current 
uncertainties in the galaxy luminosity functions. In Sec¬ 
tion 0 we show that our results at z « 10 are consistent 
with the binn ed measurements of the galaxy luminosi ty 
function ("e.g. lOesch et al.ir20l4 iBonwens et 'nil2015iill . 

The luminosity-accretion-rate relation at z ~ 6 can be 
fit with a triple power law. 


Luv = 1-5 X lO^*^ ergs“i Hz"! 



X 



M 


0.25 


3.5 Mqji 


f M 

X IH-- 

^ lOOOMeyr-i 


-0.95 

(23) 


Muv = -11.3- 2.88 log 

+ + 1000 ^ 7 ^)- 

Because of the degeneracy between the triple power law 
parameters, we choose to fix the low, middle, and high- 
M slopes at 1.15, 1.40, and 0.45, respectively based on 
Figure SI There are still degeneracies between the am¬ 
plitude and the characteristic accretion rates. The resid¬ 
uals for the fit are < 5% across the accretion rate range 
of interest. If we assume no evolution in the luminosity- 



accretion-rate relation, then Equations E51 and [Ml can be 
used as a universal EoR template to construct a fiducial 
model for the evolution of the galaxy luminosity function 
as discussed in Section 0 


6. STAR FORMATION EFFICIENCY 

To better understand the high-redshift galaxy-halo 
connection from abundance matching, we examine 
the implications of the inferred luminosity-mass and 
luminosity-accretion rate relations on the star formation 
rate and efficiency. The star formation rate Ms can be 
estimated fro m the UV luminos ity using the standard 
relation from IMadau et all (|1998f) : 


Ms = 


^uv 


8 X 10^’^ ergs^^Hz" 


Moyr” 


(25) 


assuming a Salpeter initial mass function (IMF) trun¬ 
cated at 0.1 Mq and 125 Mq. The normalization is 
uncertain and depends on the stellar IMF and forma¬ 
tion history, but we adopt this commonly used relation 
to allow a wider comparison with other work. 

We calculate the star formation efficiency in different 
ways for the two abundance matching results. In the 
case of the luminosity-accretion-rate relation, the star 
formation efficiency is defined as 


= ^ 
~ Mb’ 


where the baryonic mass accretion rate. 


Mb = §^M cx (1 + z)2 

i irn 


(26) 


(27) 


is calculated assuming the cosmic baryon fraction. In the 
case of the luminosity-mass relation, we assume no infor¬ 
mation about accretion rates and define the efficiency as 


Ms 

MbAdyn’ 


(28) 


where the baryonic mass consumption rate, 


Mb 

^dyn 



/ 32G*Phalo 


1/2 

oc (l + z)i'5. 


(29) 


depends on the redshift-dependent halo dynamical time. 
Equation [^ni can be considered the 3-dimensional analog 
of the observed Schmidt-Kennicutt relation and is often 
used in semi-analytical models and cosmological simula¬ 
tions of galaxy formation. The star formation efficiency 
is expected to range from ^ 0.01 to ~ 0 . 1 , but keep in 
mind that the normalization is uncertain. 

We impose a star formation limit based on the crite¬ 
rion that galaxies form in dark matter halos where the 
gas cools efficiently through atomic transitions. Eor a 
minimum virial temperature Tmin = 10^ K, the mini- 
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Figure 5. Left: The star formation efficiency as a function of halo mass for z ^ 6 (blue), 8 (green), and 10 (red) from abundance matching 
Luv and M. The shaded curves are illustrative of the Icr (68%) uncertainty in the currently-observed galaxy luminosity functions. Right: 
The star formation efficiency as a function of accretion rate from abundance matching L^jv and M. The results are consistent with no 
evolution for the redshift range 6 ^ z < 10. 


mum mass for a halo to host a galaxy is defined as 




\fj,mHGj V 3 


r, Phalo 


- 1/2 


« 1.6 X 10® Mq 


^halo 

200 0.132 


- 1/2 



(30) 


where the mean atomic mass is set as /i = 0.6 for a 
fully ionized intrahalo medium. For our fiducial model 
assuming a universal EoR luminosity-accretion rate re¬ 
lation, the corresponding limiting magnitude can be fit 
with the linear relation, 

Muv.sf(2) = -10.0-0.35(z-6). (31) 

Thus, we truncate all relevant curves in our figures at 
the star formation limit. 

Figure [S] shows the inferred star formation efficiencies 
from abundance matching Tuv with M and with M. The 
uncertainties in the efficiencies correspond to the those in 
the luminosity-mass and luminosity-accretion-rate rela¬ 
tions shown in Figure |4l The efficiency is not monotonic 
with mass nor accretion rate at any given redshift, but 
has a maximum value at a characteristic peak scale near 
where the galaxy luminosity function transitions from a 
power law to an exponential decline. This peak occurs at 
a characteristic mass ^ 2 x 10^^ Mq and a characteristic 
accretion rate ~ 6 x 10^ MQyr~^ at z ~ 6. 

The dependence of the star formation efficiency on 
mass and accretion rate has a physical explanation. The 
reduced efficiency at higher M and M is consistent with 
relatively inefficient atomic cooling and cold gas accre¬ 
tion in larger halos. The atomic cooling rate is relatively 


low at T > 10® K (e.g. [Sutherland fc DoDitalll993h and 
this coincides with the observed peak mass and accretion 
rate. The reduced efficiency at lower M and M is consis¬ 
tent with feedback effects from photoheating and super¬ 
nova. Smaller halos are more severely affected by feed¬ 
back because of the lower binding energy which scales as 
M®/®. It is possible that the assumed Schechter form for 
the galaxy luminosity function, in particular the bright- 
end exponential decline and the extrapolated faint-end 
power law, may not be accurate and therefore would bias 
the inferred star formation efficiency. However, there 
is support for the nonmonotonic effic iencv, which has 
been observed at lower redshifts (e.g . iGuo et al.l 120101 : 
iLeauthaud et al.ll2012t iBehroozi et ah 12013^ 

We find that the star formation efficiency inferred 

from abundance matching Tuv and M is more con¬ 
sistent with having no redshift evolution than the effi- 
cie ncy em inferred from a bundance matching Luv and 
M. IBehroozi et al.l (j2013a[l have also found that there is 
lack of evolution in the star f ormation efficiency for 
the redshift range 0 < z < 8. iFinkelstein et al.l (12015b(l 
have shown that the stellar baryon fraction increases at 
higher redshifts over the range 4 < z < 8, but this is 
consistent with both of our abundance matching results 
and does not obviously favor one particular scenario. 
Our results suggest that the star formation rate evolves 
more like (1 -I- z)^ ® rather than the commonly-assumed 
(1 -I- z)^ ® scaling from the 3-D analog of the Schmidt- 
Kennicutt relation. More precise measurements of the 
galaxy luminosity function and a more rigorous statisti¬ 
cal analysis are required to strengthen this argument. 

We note that the comparison of the exponential tails 
in the galaxy and halo abundances is very sensitive to 
the value of M* in the Schechter luminosity function. 
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In Figure [51 the star formation efficiency curves are 
shifted to higher M and M because having M* = —20.92 
at 0 ~ 10 is rather bright considering that the best- 
fit characteristic magnitude changed from -20.94 to - 
20.63 for z « 6 to 8. Recall that the extrapolation 
of the observed galaxy luminosity functions from a dif- 
ferent redshift range 5 ^ 2 : 8 yields M* « —20.22 

([Bouwens et aLll^OlSbU . With this fainter characteristic 
magnitude, we expect the peaks in the star formation 
efficiency in Figure [5] (right) to be even more aligned. 

The efhciency at z ~ 6 can be fit with a triple power 
law, 




0.15 


X l-f 


X 1-b 


M 


0.25 


3.5 MQ-yi 


M 


-0.95 


1000 MQ-yi 


-1 


(32) 


Equation [32] corresponds to the best fit for the 
luminosity-accretion-rate relation in Equation 1231 This 
universal EoR template can be used in semi-analytical 
calculations and cosmological simulations to model star 
format ion in high-red s hift g alaxies during the EoR. Re¬ 
cently, iMashian et al.l (1201511 have also applied the abun¬ 
dance matching technique to calibrate a relation between 
star formation rate and halo mass. They find a double 
power law scaling relation for M > 10^^Mq, which is 
consistent with our results when the same mass range is 
considered. 


7. GALAXY LUMINOSITY FUNCTION 



Figure 6. The galaxy luminosity function at 2 : 6 (blue), 8 

(green), and 10 (red) for our fiducial model assuming a universal 
EoR luminosity-accretion-rate relation. The binned observational 
measurements at 2 : 6 ( blue squares), 8 (green circles), and 10 

(red triangles) fromIBouwens et al.l fl2015al) and at 2 : 10 (yellow 

triangles) from lOesch et alTil20^T are shown for comparison. The 
luminosity functions match at 2 6 by construction and are in 

very good agreement at 2 8. At z ^ 10, our fiducial model is 

still consistent with the highly uncertain observations. 


Table 2 

Galaxy luminosity function parameters 


To understand how the abundance of galaxies as a 
function of luminosity and redshift is connected to the 
abundance and growth of dark matter halos, it is in¬ 
formative to relate the galaxy luminosity function to 
the halo accretion rate function and the luminosity- 
accretion-rate relation. The luminosity function in mag¬ 
nitude form can be calculated as 


z) 


dn dM 
dM dM\jY ’ 


(33) 


where dn/dM is obtained using Equation [TOl Since the 
luminosity-accretion-rate relation and the star formation 
efficiency are consistent with no evolution for the redshift 
range 6 < z < 10, it is intriguing to predict the galaxy 
luminosity function using the universal EoR template for 
M\jy{M) given by EauationIMl If more precise measure¬ 
ments of the luminosity function from upcoming obser¬ 
vations clearly deviate from this fiducial model, then it 
would point to exciting, additional astrophysics in star 
and galaxy formation. 

Figure jS] shows the predicted galaxy luminosity func¬ 
tion at z Ri 6, 8, and 10 for our fiducial model. The 
shaded prediction curves reflect the Icr (68%) uncer¬ 
tainty in the luminosity-accretion rate relation. The 
binned observa t ional m easurements at 6 z < 10 from 
iBouwens et al.l (|2015all and at z ~ 10 from lOesch et al.l 
(1201411 are shown for comparison. Our results match 


2 

<P* 

M* 

a 

6 

4.9 X 10“"' 

-20.92 

-1.88 

7 

3.6 X lO-"* 

-20.76 

-1.94 

8 

2.4 X lO--* 

-20.61 

-2.01 

9 

1.5 X lO--* 

-20.46 

-2.08 

10 

8.9 X 10“® 

-20.32 

-2.15 

11 

5.0 X 10“® 

-20.18 

-2.22 

12 

2.7 X 10-5 

-20.04 

-2.30 

13 

1.4 X 10“5 

-19.91 

-2.38 

14 

7.0 X 10-5 

-19.78 

-2.45 

15 

3.4 X 10-5 

-19.64 

-2.53 


Note. — The best-fit Schechter parameters for our fiducial 
model. The normalization </>* has units of Mpc~®. 


at z ~ 6 by construction and are in very good agree¬ 
ment at z Ri 8. At z K. 10, our model generally has 
larger amplitude than the few observational data points. 
This explains the apparent redshift dependence in the 
luminosity-accretion rate relation (Figure|T]) and the star 
formation efficiency (Figure [5]). However, our model is 
still consistent with the highly uncertain observations. 
Thus, the luminosity-accretion rate relation (Equation 
[Ml) and the star formation efficiency (Equation 15^ are 
highly consistent with no evolution an d can be used 
as uni v ersal EoR templates for the EoR. IMashian et al.l 
(|2015ll : [Mason et al.l ([201511 have also predicted similar 
evolution for the galaxy luminosity function by relating 
observed star formation rate to halo mass. 
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Table [5] shows the best-fit Schechter parameters for 
our fiducial model for the redshift range 6 < 2 < 15. 
At higher redshifts, the normalization 0* decreases dra¬ 
matically, the characteristic magnitude M* is more pos¬ 
itive (fainter), and the faint-end slope a is more nega¬ 
tive (steeper). While there are degeneracies between the 
best-fit parameters, all three must vary with redshift to 
have a good fit. Note that the differences: Alog0(2:) = 
log 0 ( 2 ) - log(/)(6), AMuy{z) = Muv(z) - Muy{6), and 
Aq!(z) = a{z) — a{6) are more general and can be ap¬ 
plied to estimate the evolution of the galaxy luminosity 
function when a different calibration at z ~ 6 is used. 
For better accuracy, we suggest repeating the abundance 
matching procedure using the new galaxy luminosity 
function and the halo accretion rate function. 

We can explain the redshift dependence of the 
Schechter parameters by considering what happens to 
the halo accretion rate function (Figure [3|) with increas¬ 
ing redshift. The overall amplitude of the halo accretion 
rate function decreases and the normalization also de¬ 
creases accordingly. The characteristic accretion rate, at 
which the halo abundance undergoes exponential decline, 
shifts to lower values. Correspondingly, the characteris¬ 
tic luminosity L* decreases and the characteristic mag¬ 
nitude M* becomes more positive (fainter). Halos at a 
given accretion rate are rarer and the slope of the halo 
accretion rate function becomes steeper. Consequently, 
the faint-end slope a becomes more negative. 

If more precise measurements of the galaxy luminos¬ 
ity function find different amplitudes at higher redshifts 
compared to our fiducial model, then it could be due to 
a number of reasons. It may simply be that the galaxy 
luminosity function at z ~ 6 used to calibrate our model 
may be overestimated or underestimated given the cur¬ 
rent measurement errors. Or perhaps it is due to addi¬ 
tional astrophysics in star and galaxy formation. Lower 
metallicity and less metal cooling could lead to a larger 
Jeans mass for molecular clouds and therefore result in a 
more top-heavy IMF. Less photoheating and less Jeans 
smoothing of the IGM at an earlier stage in reionization 
could result in more efficient gas accretion and star for¬ 
mation. This would affect the formation and abundance 
of dwarf galaxies and lead to different values for the faint- 
end slope. More accurate measurements are required to 
understand the dependence of the star formation effi¬ 
ciency on mass, accretion rate, and redshift and to shed 
light on the galaxy formation process. 

8. FORECAST FOR JWST 

Future observations with the James Webb Space Tele¬ 
scope (JWST) will provide an exciting and important 
window to the EoR. It will allow more precise measure¬ 
ments of the galaxy luminosity function, particularly at 
lower luminositie s and higher redshifts. According to 
iWindhorst et al.l (12014 and through personal communi¬ 
cation) , JWST has the sensitivity to reach AB apparent 
magnitude uiab ~ 31 with deep observations and even 
rriAB ~ 32 with ultra deep observations. They have sug¬ 
gested the strategy of observing a large number of deep 
fields and a much larger number of medium-deep surveys 
on gravitational leasing foreground targets. 

Figure [7] and Table |3] show the forecast based on our 
fiducial model for the evolution of the galaxy luminos¬ 
ity function. From z = 6toz = 15, the limiting ab- 



26 28 30 32 


Figure 7. The cumulative number of galaxies brighter than muv 
per square degree per unit redshift at 2 = 6 (blue), 10 (green), 
13 (yellow), and 15 (red) for our fiducial model. JWST has the 
sensitivity to observe unlensed galaxies at least down to M* (black 
x) at 2 : < 13 (15) with apparent magnitude limit (grey) ttiab ~ 31 
(32). 

solute magnitude Mab gets brighter by ~ 2.3 magni¬ 
tudes as the luminosity distance increases by a factor 
of ~ 2.8. The comoving number density of observable 
galaxies n{< Mab) drops by about 4 — 5 orders of mag¬ 
nitude. The number of observable galaxies per square 
degree per unit redshift dN{< MAB)/dz changes simi¬ 
larly since the differential comov ing volume dV/dz onl y 
changes by a facto r of ~ 2.2. iMashian et al.l (|20I5l) : 
iMason et al.l (|2015l) have also made similar forecasts for 
JWST using their predictions for the evolution of the 
galaxy luminosity function. 

JWST has the sensitivity to observe >11 (24) unlensed 
galaxies per square degree per unit redshift at least down 
to M* at z < 13 (14) with deep (ultra deep) observations. 
It will be able to probe some portion of the faint end at 
lower redshifts, but it is still about 8 (7) magnitudes away 
from the expected star formation limit. At z > 13 (14), 
JWST will mainly probe the exponential tail of the lu¬ 
minosity function, where the rarity of galaxies will make 
them difficult to find. With apparent magnitude limit 
ruAB ~ 32, it can actually reach M* at z ~ 15. Note 
that our forecasted counts will increase or decrease de¬ 
pending on how the actual galaxy luminosity function at 
z Ri 6 compares to the one used to calibrate our fiducial 
model. 

Our fiducial model for the evolution of the galaxy lumi¬ 
nosity function and the corresponding forecast for JWST 
are for our assumed cosmological parameters. If a differ¬ 
ent cosmological model is considered, then the results at 
z r; 6 would remain the same since this is where the 
calibration is done, but the higher redshift results would 
change. For example, consider what happens if or erg 
is increased. The ratio of the halo accretion rate function 
at redshift z > 6 relative to that at the pivot z = 6 will 
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Table 3 

Forecast of Galaxy Counts for JWST 


z 

M 31 

M 32 

n(< M 31 ) 

n(< M 32 ) 

dN(< M3i)ldz 

dN(< Mzi)ldz 

6 

- 17.9 

- 16.9 

5.1 X 10 "® 

1.3 X 10-3 

50000 

130000 

7 

- 18.3 

- 17.3 

2.4 X 10-3 

6.8 X 10-3 

21000 

60000 

8 

- 18.6 

- 17.6 

9.7 X lO-"' 

3.2 X 10-3 

7800 

26000 

9 

- 18.9 

- 17.9 

3.5 X lO-"' 

1.4 X 10-3 

2600 

9900 

10 

- 19.2 

- 18.2 

1.2 X lO-"' 

5.3 X 10-4 

760 

3500 

11 

- 19.4 

- 18.4 

3.4 X 10 “® 

1.9 X 10-4 

210 

1200 

12 

- 19.6 

- 18.6 

8.7 X 10 -® 

6.1 X 10-4 

49 

340 

13 

- 19.8 

- 18.8 

2.1 X 10 -® 

1.8 X 10-3 

11 

94 

14 

- 20.0 

- 19.0 

4.2 X 10 “’^ 

5.0 X 10-3 

2 

24 

15 

- 20.2 

- 19.2 

7.8 X 10-3 

1.2 X 10-3 

0.4 

6 


Note. — The absolute magnitude Mx corresponds to an apparent magnitude limit rriAB = Cumulative comoving number density 
n(< Mab) has units of Mpc“^ and cumulative galaxy count per angular area per unit redshift dN{< M\s)/dz has units of deg~^. 


increase. As a result, the amplitude of the galaxy lumi¬ 
nosity function and forecast for JWST at higher redshifts 
will also increase compared to our fiducial predictions. 


9. DISCUSSION 

The reionization history is determined by the evolving 
abundance of escaped ionization photons, which depends 
on the UV luminosity function, the spectral energy distri¬ 
bution (SED), and the radiation escape fraction of high- 
redshift galaxies. Only a small fraction of the ionizing 
photon budget comes from currently observable galaxies 
with Muv ^ —17 and at z < 10. In order to calculate 
the total budget, we need to extrapolate the luminos¬ 
ity function to fainter magnitudes and higher redshifts. 
With our fiducial model, we can more confidently inte¬ 
grate the luminosity function because of the well-behaved 
and physically-motivated evolution in the Schechter pa¬ 
rameters. 

Figure |8] shows our prediction for the cumulative ion¬ 
izing photon number density, 

P^max. 

n-y{> z)= n.^(z) 

J Z 

where the photon production rate density is calculated 
as 

nj(z) = / (l){M\jy,z)N^{M\jy)dMuv- (35) 

J Albright 

A galaxy with Population II st ars produces ionizing pho¬ 
tons at a rate given by (e.g. iBruzual fc Charlotl 120031 
lSchaereJl2003ll 


dz, 


(34) 




^q45.9— 0.4Muv g— 1 


1025.2 g-1 


Luv 


ergs 


-iHz-i 


(36) 


For the integration limits, we choose Zmax = 25 and 
Afbright = — 5. For Mfaint, we choose the star for¬ 

mation limit given by Equation m or the commonly- 
assumed limit Mttv = —10 (e.g. iBouwens et ^ 120121 : 
iR.obertson et al.ll2013[l for comparison. We also show the 
contribution probed by JWST with apparent magnitude 
limit muv ~ 31 — 32. 

Our fiducial model for the galaxy luminosity function 
at 6 < z < 8 is in very good agreement with current ob¬ 
servations. At z ~ 10, it generally has higher amplitude 
but shallower faint-end slope than the highly-uncertain, 



Figure 8. The cumulative ionizing photon count per hydrogen 
atom from our fiducial model. The galaxy luminosity function 
is integrated down to the star formation limit (blue) and to the 
commonly-assumed limit Muv = —10 (red dotted). Also shown is 
the contribution probed by JWST with apparent magnitude limit 
muv ~ 31 - 32 (grey). 


best-fit observations. Increasing (reducing) the ioniz¬ 
ing photon production rate at higher redshifts would 
hasten (delay) the start of reionization and lengthen 
(shorten) its duration. Subsequently, this would af¬ 
fect integrated measurements of electron scattering on 
the CMB such as patchy Thomson scattering, kinetic 
Sunyaev-Zel’dovich temperature anisotropy, and polar¬ 
ization anisotropy. However, there is still no tension with 
current EoR constraints since we poorly understand how 
the radiation escape fraction varies with galaxy luminos¬ 
ity, halo mass, and redshift. 

We impose a physically-motivated star formation limit 
based on the criterion that galaxies form in dark matter 
halos where the gas cools efficiently through atomic tran¬ 
sitions. For our fiducial model, we find Muv,SF ~ —10 
at z = 6, but it becomes more negative (brighter) 
with increasing redshift. This also reduces the ion- 




















SCORCH I 


13 


izing photon production rate at higher redshifts com¬ 
pared to calculations using the commonly-assumed limit 
Muv = —10 at all redshifts. Careful consideration is 
necessary when integrating the galaxy luminosity func¬ 
tion down to the faintest limit, especially if the faint-end 
slope is steep. For example, if we extrapolate the fit¬ 
ting form ula for the evolution of the Schechter parame¬ 
ters from iBouwens et all (|2015al) and integrate down to 
the commonly-assumed magnitude limit, then the pho¬ 
ton counts increase by more than an order of magnitude 
at high redshifts. 

The cumulative photon count per hydrogen atom 
reaches unity at z ~ 10 and therefore, reionization is 
completed below this redshift since we have yet to ac¬ 
count for the radiation escape fraction. Assuming reion¬ 
ization ended at 2 > 6 , the luminosity-weighted radia¬ 
tion escape fraction is constrained by (/esc) > 0.1. The 
lower limit is actually higher than this because addi¬ 
tional photons are required to balance recombinations 
in the clumpy IGM. Since we obs erve radiation escap e 
fractions /esc 0.01 at z 3 (e.g lShaolev et al.]l2006ll . 
this suggests that the escape fraction should vary with 
redshift. In addition, the average number of recombi¬ 
nations per hydrogen atom is < 10 because the escape 
fraction cannot exceed unity by definition. While these 
are only basic constraints, they already help to narrow 
down the parameter space of reionization. In an upcom¬ 
ing paper for the SCORCH project, we will make more 
robust constraints on the radiation escape fraction and 
hydrogen clumping factor. 

In large-scale reionization simulations with box sizes 
> 50 Mpch“^, radiation sources are often identified 
with dark matter halos from N-body simulations, with 
properties modeled using simple scaling relations. The 
source luminosity L is assumed to be proportional to 
the halo mass M thro ugh the light-to-m ass ratio, which 
can be co nstant fe.g. [ifiey et al.l [ 200611 . have mass de¬ 


pendence (|McQuinn et al 


2007), or have both mass and 


redshift dependence fe.g. iTrac fc CeDll2007D . However, 
in all three cases the luminosity is deterministic, with¬ 
out any scatter to account for the episodic nature of 
starbursting galaxies. In an upcoming paper for the 
SCORCH project, we will use a physically-motivated 
and observationally-constrained approach for modeling 
galaxy formation in large-scale reionization simulations. 
Using the luminosity-accretion rate relation, we can 
study the effects of episodic star formation on the dis¬ 
tribution and morphology of HII regions. 


10. CONCLUSIONS 

SCORCH is a new project to study the EoR and pro¬ 
vide useful theoretical tools and predictions to facilitate 
more accurate and efficient comparison between observa¬ 
tions and theory. In this first paper, we probe the con¬ 
nection between observed high-redshift galaxies and sim¬ 
ulated dark matter halos in order to better understand 
the distribution and evolution of the primary source of 
ionizing radiation. A series of 22 high-resolution N-body 
simulations shown in Table [T] is used to quantify the 
abundance of dark matter halos as a function of mass 
M, accretion rate M, and redshift z. The abundance 
matching technique is used to connect the distribution 
of observed high-redshift galaxies to the distribution of 
simulated dark matter halos. The major results are: 


1. The halo mass function dn/dM can be calcu¬ 
lated using a self-similar barrier-crossing distribu¬ 
tion function /(cr) given by Equation [H] The new 
fit is ~ 20 % more accurate at the high-mass end 
where bright galaxies are expected to reside. 

2. The distribution of mass accretion rate at any given 
mass is positively skewed. The average accretion 
rate (M), variance ^ 2 , and skewness ^3 as functions 
of mass and redshift are quantified by EauationsfTTl 
[121 anddH 

3. The halo accretion rate function dn/dM is related 
to the halo mass function through a mediating mass 
relation, and can be calculated using Equations [TUI 
and 1 141 Both halo abundance functions are similar 
in shape, with each having a low-end power law 
and a high-end exponential decline. 

4. The luminosity-accretion relation Luv(Af, z) is 
more consistent with no redshift evolution than the 
luminosity-mass relation Luy{M, z) for the red¬ 
shift range 6 < z < 10. The former is quantified 
with a universal EoR template given by Equations 

[23] and [Ml 

5. The star formation rate Ms evolves more like (1 -I- 
z)^'^ rather than the commonly-assumed (1 -I- z)^ ® 
scaling. More precise measurements of the galaxy 
luminosity function and a more rigorous statistical 
analysis are required to strengthen this argument. 

6 . The star formation efficiency e is not mono¬ 
tonic with mass nor accretion rate, but reaches a 
maximum value at a characteristic mass ^ 2 x 
10 ^^ Mq and a characteristic accretion rate ^ 
6 X 10^ M 0 yr“^ at z ~ 6 . A corresponding tem¬ 
plate for the efficiency as a function of accretion 
rate is given by Equation 15^ 

7. The faintest magnitude corresponding to the star 
formation limit is Muv.SF ~ —10 at z = 6, but 
it becomes more negative (brighter) with increas¬ 
ing redshift. The redshift dependence is given by 
Equation [STJ 

8 . The fiducial model for the galaxy luminosity func¬ 
tion is constructed from the halo accretion rate 
function and the luminosity-accretion-rate relation 
using Equation |33l The evolution of the Schechter 
parameters are shown in Table [U </)* decreases, M* 
is more positive (fainter), and a is more negative 
(steeper) at higher redshifts. 

9. Forecast galaxy counts for JWST are shown in Ta¬ 
ble |3l JWST has the sensitivity to observe >11 

(24) unlensed galaxies per square degree per unit 
redshift at least down to M* at z < 13 (14) with 
deep (ultra deep) observations. It will also be able 
to probe some portion of the faint end at lower 
redshifts. 

The numerical fits for the abundance matching re¬ 
sults: luminosity-accretion-rate relation, star formation 
efficiency, Schechter luminosity function parameters, and 
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JWST forecasts are based on the current galaxy luminos¬ 
ity function at z « 6. When updated observations are 
available, we suggest repeating the abundance matching 
procedure using the new galaxy luminosity function and 
the halo accretion rate function. 
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